"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBRead.py 990 2012-02-15 13:52:11Z bgoli $)

"""

import os, time
import numpy
##  cDir = os.path.dirname(os.path.abspath(os.sys.argv[0]))
import CBXML
from CBVersion import __DEBUG__, __version__

##  try:
    ##  import psyco
    ##  psyco.full()
##  except:
    ##  pass


_HAVE_SYMPY_ = None
try:
    import sympy
    _HAVE_SYMPY_ = True
except ImportError:
    print 'Rational IO not available'
    _HAVE_SYMPY_ = False
_HAVE_SYMPY_ = None
try:
    import h5py
    _HAVE_HD5_ = True
except ImportError:
    _HAVE_HD5_ = False
    

def readSBML3(fname, work_dir=None):
    """
    Read in an SBML L3 FBA file where:

     - *fname* is the filename
     - *work_dir* is the working directory

    Note this functions is a backwards compatible wrapper for CBXML.readSBML3
    """
    print '\n*****\nThis function is being deprecated, please use readSBML2FBA() instead!\n*****\n'
    return CBXML.sbml_readSBML2FBA(fname, work_dir, return_sbml_model=False)

def readSBML2FBA(fname, work_dir=None, return_sbml_model=False):
    """
    Read in an SBML Level 2 file with FBA annotation where:

     - *fname* is the filename
     - *work_dir* is the working directory
     - *return_sbml_model* [default=False] return a a (cbm_mod, sbml_mod) pair

    """    
    return CBXML.sbml_readSBML2FBA(fname, work_dir, return_sbml_model)

def readLPtoList(fname, work_dir):
    NEW = False
    TYPE = None
    Object = []
    Constr =[]
    Bounds = []
    F = file(os.path.join(work_dir, fname), 'r')
    for l in F:
        if l == '' or l[:2] == '\\\\' or l == '\n' or l.strip() == 'END':
            print 'skipping'
            if __DEBUG__: print l
        else:
            L = l.strip()
            if L == 'Maximize':
                TYPE = 'ObjFunc'
                NEW = True
            elif L == 'Subject To':
                TYPE = 'Constr'
                NEW = True
            elif L == 'Bounds':
                TYPE = 'Bounds'
                NEW = True
            else:
                NEW = False

            if TYPE == 'ObjFunc' and not NEW:
                Object.append(L)
            elif TYPE == 'Constr' and not NEW:
                Constr.append(L)
            elif TYPE == 'Bounds' and not NEW:
                Bounds.append(L)

    F.close()
    if __DEBUG__:
        print 'ObjectiveLines'
        print Object
        print 'ConstraintLines'
        print Constr
        print 'BoundsLines'
        print Bounds
    return Object, Constr, Bounds

def readSK_FVA(filename):
    """
    Read Stevens FVA results (opt.fva) file and return a list of dictionaries

    """
    assert os.path.exists(filename), '\nGive me a break!\n'
    ##  name = []
    vari = []
    F = file(filename, 'r')
    for l in F:
        L = l.split(':')
        Jn = L[0].strip()
        Vmin = L[1].strip()
        V = L[2].strip()
        V = V.split('--')
        Vmax = V[0].strip()
        Vstat = V[1].strip()
        ##  name.append((Jn, Vstat))
        ##  vari.append((Vmin, Vmax))
        vari.append({'name' : Jn,
                     'min' : Vmin,
                     'max' : Vmax,
                     'status' : Vstat
                    })
        if __DEBUG__: print Jn, Vmin, Vmax, Vstat
    return vari

##  def readSK_vertexOld(fname, bigfile=False):
    ##  """
    ##  Reads in Stevens vertex analysis file and returns:

     ##  - a list of vertex vectors
     ##  - a list of ray vectors
     ##  - the basis of the lineality space as a list of vectors

    ##  all vectors in terms of the column space of N

    ##  """

    ##  assert _HAVE_SYMPY_, 'Install Sympy for rational IO support'

    ##  assert os.path.exists(fname), 'Uhm exqueeze me ...'
    ##  SK_vert_file = file(fname, 'r')
    ##  VertOut = []
    ##  LinOut = []
    ##  RayOut = []
    ##  if bigfile:
        ##  VertTmp = gzip.open('_vtx_.tmp.gz','wb', compresslevel=3)
        ##  LinTmp = gzip.open('_lin_.tmp.gz','wb', compresslevel=3)
        ##  RayTmp = gzip.open('_ray_.tmp.gz','wb', compresslevel=3)
    ##  GOvert = False
    ##  GOray = False
    ##  GOlin = False
    ##  lcntr = 0
    ##  lcntrtmp = 0
    ##  for l in SK_vert_file:
        ##  lcntr += 1
        ##  if lcntr == 1000:
            ##  print 'Processing vertex: %s' % (lcntr + lcntrtmp)
            ##  lcntrtmp += lcntr
            ##  lcntr = 0
        ##  if '* Lineality basis ' in l:
            ##  GOvert = False
            ##  GOray = False
            ##  GOlin = True
        ##  if '* Rays ' in l:
            ##  GOvert = False
            ##  GOray = True
            ##  GOlin = False
        ##  if '* Vertices ' in l:
            ##  GOvert = True
            ##  GOray = False
            ##  GOlin = False
        ##  if l[:2] != '* ':
            ##  L = l.split()
            ##  rowL = []
            ##  for c in L:
                ##  c = c.strip()
                ##  if c == '0':
                    ##  rnum = '0'
                    ##  rowL.append(0.0)
                ##  else:
                    ##  rnum = sympy.Rational('%s' % c)
                    ##  rowL.append(rnum.evalf())
                ##  del rnum
            ##  del L
            ##  if GOlin:
                ##  if bigfile:
                    ##  rowEnd = len(rowL)
                    ##  cntr = 0
                    ##  for e in rowL:
                        ##  cntr += 1
                        ##  if e == 0.0:
                            ##  LinTmp.write('0.0')
                        ##  else:
                            ##  LinTmp.write('%.14f' % e)
                        ##  if cntr == rowEnd:
                            ##  LinTmp.write('\n')
                        ##  else:
                            ##  LinTmp.write(',')
                ##  else:
                    ##  LinOut.append(rowL)
            ##  elif GOray:
                ##  if bigfile:
                    ##  rowEnd = len(rowL)
                    ##  cntr = 0
                    ##  for e in rowL:
                        ##  cntr += 1
                        ##  if e == 0.0:
                            ##  RayTmp.write('0.0')
                        ##  else:
                            ##  RayTmp.write('%.14f' % e)
                        ##  if cntr == rowEnd:
                            ##  RayTmp.write('\n')
                        ##  else:
                            ##  RayTmp.write(',')
                ##  else:
                    ##  RayOut.append(rowL)
            ##  elif GOvert:
                ##  if bigfile:
                    ##  rowEnd = len(rowL)
                    ##  cntr = 0
                    ##  for e in rowL:
                        ##  cntr += 1
                        ##  if e == 0.0:
                            ##  VertTmp.write('0.0')
                        ##  else:
                            ##  VertTmp.write('%.14f' % e)
                        ##  if cntr == rowEnd:
                            ##  VertTmp.write('\n')
                        ##  else:
                            ##  VertTmp.write(',')
                ##  else:
                    ##  VertOut.append(rowL)
            ##  del rowL
    ##  print '\nProcessed %s vertices.\n' % (lcntr + lcntrtmp)
    ##  SK_vert_file.close()
    ##  if bigfile:
        ##  VertTmp.close()
        ##  RayTmp.close()
        ##  LinTmp.close()
        ##  VertTmp = gzip.open('_vtx_.tmp.gz','rb')
        ##  LinTmp = gzip.open('_lin_.tmp.gz','rb')
        ##  RayTmp = gzip.open('_ray_.tmp.gz','rb')        
        ##  return VertTmp, RayTmp, LinTmp
    ##  else:
        ##  print 'Lineality basis: %s' % len(LinOut)
        ##  print 'Number of rays: %s' % len(RayOut)
        ##  print 'Number of vertices: %s' % len(VertOut)
        ##  return VertOut, RayOut, LinOut

def readSK_vertexOld(fname, bigfile=False, fast_rational=False, nformat='%.14f', compresslevel=3):
    """
    Reads in Stevens vertex analysis file and returns, even more optimized for large datasets than the original.

     - a list of vertex vectors
     - a list of ray vectors
     - the basis of the lineality space as a list of vectors

    all vectors in terms of the column space of N

    """

    if fast_rational:
        pass
    else:
        assert _HAVE_SYMPY_, 'Install Sympy for rational IO support'

    assert os.path.exists(fname), 'Uhm exqueeze me ...'
    
    print '\n**********\nreadSK_vertex options are:\n'
    print 'bigfile: %s' % bigfile
    print 'fast_rational: %s' % fast_rational
    print 'nformat: %s' % (nformat % 0.12345678901234567890)
    print '**********\n'
    
    SK_vert_file = file(fname, 'r')
    VertOut = []
    LinOut = []
    RayOut = []
    if bigfile:
        VertTmp = gzip.open('_vtx_.tmp.gz','wb', compresslevel=compresslevel)
        LinTmp = gzip.open('_lin_.tmp.gz','wb', compresslevel=compresslevel)
        RayTmp = gzip.open('_ray_.tmp.gz','wb', compresslevel=compresslevel)
    GOvert = False
    GOray = False
    GOlin = False
    lcntr = 0
    lcntrtmp = 0
    TZero = time.time()
    for l in SK_vert_file:
        lcntr += 1
        if lcntr == 1000:
            print 'Processing vertex: %s (%s min)' % (lcntr + lcntrtmp, round((time.time()-TZero)/60.0,1))
            lcntrtmp += lcntr
            lcntr = 0
        if '* Lineality basis ' in l:
            GOvert = False
            GOray = False
            GOlin = True
        if '* Rays ' in l:
            GOvert = False
            GOray = True
            GOlin = False
        if '* Vertices ' in l:
            GOvert = True
            GOray = False
            GOlin = False
        if l[:2] != '* ':
            L = l.split()
            rowL = []
            for c in L:
                rnum = None
                c = c.strip()
                if c == '0':
                    rnum = '0'
                    rowL.append(0.0)
                else:
                    if not fast_rational:
                        rnum = sympy.Rational('%s' % c)
                        rowL.append(rnum.evalf())
                        ##  print c, rnum.evalf()
                    else:
                        rnum = c.split('/')
                        if len(rnum) == 1:
                            rowL.append(float(rnum[0]))
                            ##  print c, float(rnum[0])
                        else:
                            rowL.append(float(rnum[0])/float(rnum[1]))
                            ##  print c, float(rnum[0].strip())/float(rnum[1].strip())
                del rnum
            del L
            rowL = tuple(rowL)
            if GOlin:
                if bigfile:
                    rowEnd = len(rowL)
                    cntr = 0
                    for e in rowL:
                        cntr += 1
                        if e == 0.0:
                            LinTmp.write('0.0')
                        else:
                            LinTmp.write(nformat % e)
                        if cntr == rowEnd:
                            LinTmp.write('\n')
                        else:
                            LinTmp.write(',')
                else:
                    LinOut.append(rowL)
            elif GOray:
                if bigfile:
                    rowEnd = len(rowL)
                    cntr = 0
                    for e in rowL:
                        cntr += 1
                        if e == 0.0:
                            RayTmp.write('0.0')
                        else:
                            RayTmp.write(nformat % e)
                        if cntr == rowEnd:
                            RayTmp.write('\n')
                        else:
                            RayTmp.write(',')
                else:
                    RayOut.append(rowL)
            elif GOvert:
                if bigfile:
                    rowEnd = len(rowL)
                    cntr = 0
                    for e in rowL:
                        cntr += 1
                        if e == 0.0:
                            VertTmp.write('0.0')
                        else:
                            VertTmp.write(nformat % e)
                        if cntr == rowEnd:
                            VertTmp.write('\n')
                        else:
                            VertTmp.write(',')
                else:
                    VertOut.append(rowL)
            del rowL
    print '\nProcessed %s vertices.\n' % (lcntr + lcntrtmp)
    SK_vert_file.close()
    if bigfile:
        VertTmp.close()
        RayTmp.close()
        LinTmp.close()
        try:
            VertTmp = gzip.open('_vtx_.tmp.gz','rb')
            LinTmp = gzip.open('_lin_.tmp.gz','rb')
            RayTmp = gzip.open('_ray_.tmp.gz','rb')        
            return VertTmp, RayTmp, LinTmp
        except Exception, ex:
            print ex
            print '\nReturning file names:'
            return '_vtx_.tmp.gz', '_lin_.tmp.gz', '_ray_.tmp.gz'
    else:
        print 'Lineality basis: %s' % len(LinOut)
        print 'Number of rays: %s' % len(RayOut)
        print 'Number of vertices: %s' % len(VertOut)
        return VertOut, RayOut, LinOut
    
def readSK_vertex(fname, bigfile=True, fast_rational=False, nformat='%.14f', compression=None, hdf5file=None):
    """
    Reads in Stevens vertex analysis file: 
    
     - *fname* the input filename (.all file that results from Stevens pipeline)
     - *bigfile* [default=True] this option is now always true and is left in for backwards compatability
     - *fast_rational* [default=False] by default off and uses SymPy for rational-->float conversion, when on uses float decomposition with a slight (2th decimal) decrease in accuracy
     - *nformat* [default='%.14f'] the number format used in output files
     - *compression* [default=None] compression to be used in hdf5 files can be one of [None, 'lzf', 'gz?', 'szip']
     - *hdf5file* [default=None] if None then generic filename '_vtx_.tmp.hdf5' is uses otherwise <hdf5file>.hdf5
    
    and returns an hdf5 *filename* of the results with a single group named **data** which countains datasets
     
     - vertices
     - rays
     - lin
     
    where all vectors are in terms of the column space of N.

    """
    
    bigfile=True

    if not fast_rational:
        assert _HAVE_SYMPY_, 'Install Sympy for rational IO support'
    if bigfile:
        assert _HAVE_HD5_, 'Install h5py for large dataset support'

    assert os.path.exists(fname), 'Uhm exqueeze me ...'
    
    print '\n**********\nreadSK_vertex options are:\n'
    print 'bigfile: %s' % bigfile
    print 'fast_rational: %s' % fast_rational
    print 'nformat: %s' % (nformat % 0.12345678901234567890)
    print '**********\n'
    
    SK_vert_file = file(fname, 'r')
    VCNTR = 0
    LCNTR = 0
    RCNTR = 0
    CCNTR = 0
    for l in SK_vert_file:
        if '* Vertices (' in l:
            VCNTR = long(l.replace('* Vertices (','').replace('vectors):','').strip())
            break
    SK_vert_file.seek(0)
    for l in SK_vert_file:
        if '* Lineality basis (' in l:
            LCNTR = long(l.replace('* Lineality basis (','').replace('vectors):','').strip())
            break
    SK_vert_file.seek(0)            
    for l in SK_vert_file:
        if '* Rays (' in l:
            RCNTR = long(l.replace('* Rays (','').replace('vectors):','').strip())
            break
    SK_vert_file.seek(0)
    for l in SK_vert_file:
        if l[:2] != '* ':
            L = l.split()
            CCNTR = len(L)
            del L
            break
    SK_vert_file.seek(0)
    print RCNTR, LCNTR, VCNTR,  CCNTR        
    VertOut = []
    LinOut = []
    RayOut = []
    outFileName = '_vtx_.tmp.hdf5'
    if bigfile:
        if hdf5file != None:
            outFileName = hdf5file+'.'+str(compression)+'.hdf5'
        HD5out = h5py.File(outFileName,'w')
        if 'data' in HD5out.keys():
            del HD5out['data']
        Dgrp = HD5out.create_group('data')
        VertTmp = HD5out['data'].create_dataset('vertices', (VCNTR, CCNTR), dtype=numpy.double, compression=compression)
        LinTmp = HD5out['data'].create_dataset('lin', (LCNTR, CCNTR), dtype=numpy.double, compression=compression)
        RayTmp = HD5out['data'].create_dataset('rays', (RCNTR, CCNTR), dtype=numpy.double, compression=compression)
    GOvert = False
    GOray = False
    GOlin = False
    lcntr = 0
    lcntrtmp = 0
    TZero = time.time()
    vert_count = 0
    lin_count = 0
    ray_count = 0
    ##  data_row = numpy.zeros((1,CCNTR), 'd')
    print '\nStarting vertex mapping at %s\n' % time.strftime('%H:%M:%S')
    for l in SK_vert_file:
        lcntr += 1
        if lcntr == 1000:
            print 'Processing vertex: %s (%3.2f at %s min)' % (lcntr + lcntrtmp, (float(lcntr + lcntrtmp)/float(VCNTR)*100.0), round((time.time()-TZero)/60.0,1))
            lcntrtmp += lcntr
            lcntr = 0
        if '* Lineality basis ' in l:
            GOvert = False
            GOray = False
            GOlin = True
        if '* Rays ' in l:
            GOvert = False
            GOray = True
            GOlin = False
        if '* Vertices ' in l:
            GOvert = True
            GOray = False
            GOlin = False
        if l[:2] != '* ':
            ##  data_row[0,:] = 0.0
            L = l.split()
            if GOlin:
                LinTmp[lin_count] = 0.0
                for c in xrange(CCNTR):
                    rnum = None
                    val = L[c].strip()
                    if val == '0':
                        ##  LinTmp[lin_count,c] = 0.0
                        pass
                    else:
                        if not fast_rational:
                            rnum = sympy.Rational('%s' % val)
                            LinTmp[lin_count,c] = rnum.evalf()
                        else:
                            rnum = val.split('/')
                            if len(rnum) == 1:
                                LinTmp[lin_count,c] = float(rnum[0])
                            else:
                                LinTmp[lin_count,c] = float(rnum[0])/float(rnum[1])
                lin_count += 1                     
            elif GOray:
                RayTmp[ray_count] = 0.0
                for c in xrange(CCNTR):
                    rnum = None
                    val = L[c].strip()
                    if val == '0':
                        ##  RayTmp[ray_count,c] = 0.0
                        pass
                    else:
                        if not fast_rational:
                            rnum = sympy.Rational('%s' % val)
                            RayTmp[ray_count,c] = rnum.evalf()
                        else:
                            rnum = val.split('/')
                            if len(rnum) == 1:
                                RayTmp[ray_count,c] = float(rnum[0])
                            else:
                                RayTmp[ray_count,c] = float(rnum[0])/float(rnum[1])
                ray_count += 1  
            elif GOvert:
                VertTmp[vert_count] = 0.0
                for c in xrange(CCNTR):
                    ##  print 'lin_count',lin_count,c
                    ##  print 'ray_count',ray_count,c
                    ##  print 'vert_count',vert_count,c
                    rnum = None
                    val = L[c].strip()
                    if val == '0':
                        ##  VertTmp[vert_count,c] = 0.0
                        pass
                    else:
                        if fast_rational:
                            rnum = val.split('/')
                            if len(rnum) == 1:
                                VertTmp[vert_count,c] = float(rnum[0])
                            else:
                                VertTmp[vert_count,c] = float(rnum[0])/float(rnum[1])                        
                        else:
                            rnum = sympy.Rational('%s' % val)
                            VertTmp[vert_count,c] = rnum.evalf()

                vert_count += 1  
                        
    print '\nProcessed %s vectors.\n' % (lcntr + lcntrtmp)
    SK_vert_file.close()
    if bigfile:
        HD5out.close()
        ##  HD5out = h5py.File(outFileName,'r')
        return outFileName
    else:
        print 'Lineality basis: %s' % len(LinOut)
        print 'Number of rays: %s' % len(RayOut)
        print 'Number of vertices: %s' % len(VertOut)
        return VertOut, RayOut, LinOut
    
if __name__ == '__main__':
    fm = readSBML3('branch.xml', work_dir='/home/bgoli/mypysces/pysces_mech/examples')
    ##  fm = readSBML3('branch.xml', work_dir='d:\\mypysces\\pysces_mech2\\examples')
    ##  fm = readSBML3('iJR904.xml.converted.xml', work_dir='/home/bgoli/mypysces/pysces_mech/examples')
    SD = fm.species
    RD = fm.reactions
    CD = fm.flux_bounds
    OD = fm.objectives

    print '\nSpecies'
    print [s.getPid() for s in SD]
    print [s.is_boundary for s in SD]
    print [s.name for s in SD]
    print [s.value for s in SD]
    print '\nReactions'
    print [r.getPid() for r in RD]
    print [r.name for r in RD]
    print [r.reversible for r in RD]
    print [r.getReagentRefs() for r in RD]
    print [r.getStoichiometry() for r in RD]
    print '\nConstraints'
    print [c.getPid() for c in CD]
    print [c.reaction for c in CD]
    print [c.operation for c in CD]
    print [c.value for c in CD]
    print '\nObjective functions'
    print [o.getPid() for o in OD]
    print [o.operation for o in OD]
    print '\t', [o.fluxObjectives[0].getPid() for o in OD]
    print '\t', [o.fluxObjectives[0].reaction for o in OD]
    print '\t', [o.fluxObjectives[0].coefficient for o in OD]
    print '\nIds from model'
    print 'fm.getPid()', fm.getPid()
    print 'fm.N', fm.N
    print fm.getSpeciesIds()
    print fm.getReactionIds()
    print fm.getFluxBoundIds()
    print fm.activeObjIdx
    print fm.setActiveObjective('J3'), fm.activeObjIdx
    print fm.setActiveObjective('J4'), fm.activeObjIdx
    print fm.getObjectiveIds()



